Anomalous Curie response of impurities 
in quantum-critical spin-1/2 Heisenberg antiferromagnets 
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We consider a magnetic impurity in two different S = 1/2 Heisenberg bilayer antiferromagnets at 
their respective critical inter-layer couplings separating Neel and disordered ground states. We cal- 
culate the impurity susceptibility using a quantum Monte Carlo method. With intra-layer couplings 
in only one of the layers (Kondo lattice), we observe an anomalous Curie constant C* , as pre- 
dicted on the basis of field-theoretical work [S. Sachdev et al., Science 286, 2479 (1999)]. The value 
C* = 0.262 ± 0.002 is larger than the normal Curie constant C = S(S + 1) /3. Our low-temperature 
results for a symmetric bilayer are consistent with a universal C* . 

PACS numbers: 75.10.Jm, 75.10.Nr, 75.40.Cx, 75.40.Mg 
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Investigating effects of various types of impurities is 
a promising strategy for studying the electronic struc- 
ture of strongly correlated systems [HQ. In very small 
concentrations, detectable changes in response functions 
cuaused by the impurities probe the inherent bulk prop- 
erties of the host material and can give importants in- 
sights complementing information derived from direct 
studies of the bulk. This strategy was recently examined 
in a field-theoretical study [1, 13] of a nearly quantum- 
critical two-dimensional (2D) antiferromagnetic host sys- 
tem with a single localized spin-S* impurity. Detailed pre- 
dictions for the magnetic response of the impurity were 
made in the different finite-temperature scaling regimes 
associated with the quantum phase transition occuring in 
the host (as a function of some coupling constant g); from 
a quantum disordered spin-gapped paramagnet (g > g c ) 
to a gapless antiferromagnet (g < g c ). The T — ► 
asymptotic impurity susceptibility — defined as the dif- 
ference between the susceptibilities of systems with and 
without impurity — was predicted to have a Curie form, 
Ximp ~ > C /T, with the Curie constant C taking different 
values depending on the coupling 

C = S 2 /3, g < g c , (1) 

C = C* = S(S+ 1)/3, S^S, g = g c , (2) 
C = S(S+l)/3, g>g c , (3) 

where we have set % = ks = 1- The classical-like re- 
sponse for g < g c , which can be understood as being due 
to the impurity spin aligning with a large cluster of corre- 
lated spins (exponentially divergent, thus with essentially 
classical dynamics as T — ► 0), was confirmed in a recent 
numerical study Q , however only up to a logarithmic cor- 
rection [formally resulting in a log-divergent C(T)\. The 
log-correction was subsequently derived quantitatively in 
a different field theoretical formulation [4j (and qualita- 
tively also using spin wave theory 0]). It is given in 
terms of known ground-state constants of the host and 
is in complete agreement with the numerical results Q- 



In the paramagnetic phase the host response is exponen- 
tially small due to the spin gap, and the usual Curie pref- 
actor S(S + l)/3 is due solely to the localized impurity 
moment. 

The most remarkable prediction of Ref. [3j is the 
anomalous Curie constant ^ in the quantum critical 
regime, i.e., for T < p s , A, where p s is the spin stiffness 
in the Neel phase and the spin gap in the paramagnetic 
phase H (both of which vanish continuously at g c ). It 
was suggested that S 2 /3 < C* < S(S + l)/3 and that 
this could be interpreted in terms of a fractional impu- 
rity spin; S 7^ S. However, the anomalous Curie constant 
was challenged by Sushkov Q , on the basis of a Green's 
function theory giving S = S also in the quantum crit- 
ical regime [lfj. The only apparent way to settle this 
issue is by explicit unbiased numerical computation of 
the impurity susceptibility of a quantum-critical system. 
To our knowledge, the only attempt so far is by Troyer, 
who carried out a quantum Monte Carlo (QMC) study 
of the 5=1/2 Heisenberg bilayer with a vacancy (effec- 
tively corresponding to an S = 1/2 impurity) 11]. This 
calculation yelded S — S within statistical errors of a 
few percent, and thus, based on this study, the anomaly 
is either not present or is very small. 

In this Letter we present results of large-scale numeri- 
cal efforts to resolve the controversy over the existence of 
an effectively fractional impurity spin at the quantum- 
critical point. We also obtain further insights into the 
temperature dependence of the impurity susceptibility 
(corrections to the Curie form) . We will present evidence 
of an anomalous Curie constant in the case of an S = 1/2 
impurity. Our result is C* = 0.262 ± 0.002, which falls 
outside the range S 2 /'3 < C* < S(S + l)/3. It should 
be noted, however, that the e-expansion presented for C* 
in Ref. [3j was not evaluated explicitly — the actual value 
was only conjectured to fall within the above range. Thus 
there is no contradiction. For an S = 1 impurity we can- 
not detect an anomaly within statistical errors, indicating 
that C* approaches S(S + l)/3 with increasing S. 
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We use an efficient stochastic series expansion (SSE) 
QMC technique [3] to compute the susceptibility of a 
single static 5 = 1/2 or 5 = 1 impurity in a quantum 
critical host. We consider two different host systems — a 
bilayer with intra- and inter-plane couplings J and Jj_, 
respectively, and an "incomplete bilayer" in which the 
J-coupling is present only in one of the layers (a Kondo 
lattice [121]). The lattices and the ways in which we in- 
troduce the impurities in them are illustrated in Fig. Q] 
The Heisenberg hamiltonian for both host systems can 
be written as 

H = J S ' Sii • Si j + A J ^ ' S2,i: • 

(hi) (hi) 

+ J±J2 S hi-S2,i, (4) 

i 

where denotes a pair of nearest neighbors on a pe- 
riodic square lattice with L x L x 2 sites and S 0) j is the 
usual spin-i operator at site i on layer a — 1,2. We define 
the tuning parameter as g — Jx/J, and J, Jj_ > 0. The 
symmetric and incomplete bilayers correspond to A = 1 
and A = 0, respectively. The quantum-critical points for 
these models have recently been extracted to high preci- 
sion; g c = 2.5220(3) and g c = 1.3888(1) for the full and 
icomple bilayers, respectively [HI]. The ground state is 
Neel ordered for g < g c and disordered (spin-gapped) for 
g > g c . As shown in Fig. [TJ by removing a single spin 
from the incomplete (a) and symmetric (b) bilayers we 
effectively introduce 5 = 1/2 impurities, as the remain- 
ing spin in the opposite layer is unpaired. An 5 = 1 
impurity (c) is obtained by making the five bonds to a 
given spin ferromagnetic (we here keep the magnitudes 
of these interactions unchanged) . 

Our investigations proceed much in the same way as 
Ref. [HI], but we have pushed the simulations to higher 
precision, larger lattices, and lower temperatures than 
previously. The calculations are very demanding due to 
the fact that the impurity susceptibility is the difference 
between two extensive quantities; the total susceptibili- 
ties with and without impurity. They are defined as 



(5) 



with Nq — 2L 2 when there is no impurity and iVi = 
11? - 1 or Ni = 2L 2 for the 5 = 1/2 and 5 = 1 impu- 
rities, respectively. The impurity susceptibility Ximp = 
Xi ~ Xo- Very large lattices are required in order to elim- 
inate finite-size effects at low temperatures — we here re- 
port results for up to 256 x 256 x 2 spins. We achieved rel- 
ative statistical errors for xo and Yi as small as « 10 ~ 6 . 
The use of improved estimators [15[ is crucial, but still 
very long simulations are required. The calculations re- 
ported here required approximately 5 x 10 5 Pentium III 
(w 1GHz) CPU hours. 





FIG. 1: (Color online) Spin impurity models. The circles 
represent S =1/2 spins with nearest-neighbor interactions J 
(in-plane) and J± (inter- plane). A vacancy in the incomplete 
bilayer (a) or full bilayer (b) host constitutes, effetively, an 
5 = 1/2 impurity. In (c), changing the signs (J — ■> — J and 
J± — > — J±) of the interactions (the five dotted bonds) of one 
of the spins results in an 5 = 1 impurity. 



Finite-size effects are investigated by considering sys- 
tems, at their quantum-critical couplings, of increasing 
size L = 4, 8, 16, . . . at each temperature T, up to L suf- 
ficiently large for ramaining finite-zise corrections to be 
negligible. Finally, size-converged results are extrapo- 
lated to zero temperature to obtain the quantum-critical 
Curie constant defined in Eq. ((2]); 



C* = lim T\ 



imp ' 



(0) 



The Curie form can be expected to apply strictly only 
in the limit T — > 0, and in practice we have to analyze 
corrections to extract C* . In Ref. ll|, finite-size scaled 



data for Txf mp were linearly extrapolated to zero tem- 
perature. We have also found an asymptotic linear cor- 
rection in all cases studied, but we disagree with Ref. 1 1 1 
in regards to the temperature at which the linear form 
is valid. In the case of a vacancy in the incomplete bi- 
layer, our results extrapolate clearly to a value for C* 
a few percent larger than 5(5 + l)/3, thus confirming 
an anomalous Curie constant (falling outside the conjec- 
tured range, however). For the symmetric bilayer, which 
is the host system previously studied by Troyer [11(, the 
linearity of the correction to C* is established only at 
the lowest temperatures we have reached, and we cannot 
reliably extrapolate to obtain a direct independent con- 
firmation of a universal (for given S) Q anomalous Curie 
constant from these results. However, the available low-T 
data points are consistent with the results for the incom- 
plete bilayer. For an 5 = 1 impurity we also find a linear 
correction to the Curie form but the extrapolated C* is 
very close to the normal value. 

We now present the data underlying our conclusions 
summarized above. In Fig. [U[a) we show the tempera- 
ture dependence of Txf mp for the incomplete bilayer, for 
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FIG. 2: (Color online) Impurity susceptibility multiplied by 
the temperature for a vacancy in the quantum-critical incom- 
plete bilayer. (a) shows all our finite-size data on a log-lin 
scale and (b) shows size converged results on a lin-lin scale. 
Examples of finite-size effects are shown in the inset. The 
solid line in (b) and the curve in (a) show a linear in T fit to 
the low-T data. 



all the system sizes we have considered and using a log- 
lin scale. Note first that as T j J — > oo the spins become 
independent moments exhibiting normal Curie behavior, 
and thus T%f mp — > —1/4, due to there being one less 
spin in the system with a vacancy than in the intact sys- 
tem. We here focus on lower T. In the limit T — ► 0, 
we have to obtain Tytf — ► +1/4 for any L (seen in 
the figure only for L = 4), due to the 5=1/2 ground 
state of the even-T systems with a vacancy. The tem- 
perature at which this can be observed is expected to 
scale as 1/L, reflecting the low-energy level spacing of a 
uantum-critical system with dynamic exponent z = 1 
Interestingly, for the larger lattices the approach to 
the limiting T = value is preceded by a miniumum 
and a maximum. The finite-size behavior for a fixed low 
T is thus also non-monotonic, as shown in the inset of 
Fig. QJb). We consider Txf mp size-converged when re- 
sults for system sizes L and L/2 agree to within statistical 
errors. We have checked this carefully for small systems 
and do not see any indications of non-monotonicities be- 
yond a single minimum. Using this criterion we have size- 
converged data for temperatures down to T/J = 1/64. 
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FIG. 3: (Color online) T times the impurity susceptibility 
for a vacancy in the quantum-critical symmetric bilayer. The 
solid line shows the linear fit of Fig. &b). The dashed (red) 
line is the linear fit of Ref. (llj. The inset shows examples of 
the size dependence. 



At T/J = 1/(64^2), the infinite -size result is probably 
marginally higher than our L = 256 result. In Fig. Hfb) 
we show the size converged results [and the almost con- 
verg cd T/J = l/(64-y/2) data] on a lin-lin scale. A lin- 
ear T dependence can be seen below T/J s» 0.12 and 
an extrapolation gives C* — 0.262(2). Note that even 
without an extrapolation it is clear that C* exceeds the 
normal Curie constant C = 1/4 because the last three 
data points in Fig. \V[b) are all above this value. 

To check for potential systematic errors, we have car- 
ried out extensive tests of our SSE codes, comparing re- 
sults based on very long runs for 4x4 systems with exact 
diagonalization data and using two independently writ- 
ten programs with different random number generators 
for large lattices. Recent calculations using a completely 
different method have also confirmed the absence of sys- 
tematic errors in SSE calculations (l6| . Finally we have 
also verified that small changes in the coupling g, of the 
order of the statistical errors of g c d3|, do not appre- 
ciably change Ximp in the range of tempertures studied 
here. We therefore consider our conclusion of an anoma- 
lous Curie constant to be beyond reasonable doubt. 

In Fig. [3] we show data for the vacancy in a complete 
bilayer; the same system that was previously studied by 
Troyer [ll| . The computational effort involved in simu- 
lating this system (for given T, L) is roughly twice that 
for the incomplete bilayer, because of the higher critical 
coupling and larger number of interactions. In addition, 
the size-convergence is slower. We have therefore not 
been able to go to as low temperatures as in the pre- 
ceding case. We also find much larger corrections to the 
linear low-T behavior, and as a consequence we cannot 
carry out a reliable T —> extrapolation. However, as 
shown in Fig. [3l the line fit to the low-T results for the 
incomplete bilayer in Fig. [2] describes reasonably well also 
the low-T results for the symmetric bilayer. The L = 256 
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FIG. 4: (Color online) Impurity susceptibility of an S = 1 
vacancy in the quantum-critical symmetric bilayer. The inset 
shows size converged results and a linear fit to the low-T data. 



result at the lowest temperature is likely not completely 
size converged. The L = 256 point at the next-lowest 
T is also below the line, but the deviation is only two 
error bars and also here there may be some remaining 
finite-size effects. Moreover, the slope is most likely not 
universal and should thus depend on the low-energy pa- 
rameters of the bulk system. The line shown may therefor 
not be the correct one for this system although C* should 
be universal (for given S) 0, |j| ■ Thus we regard these 
symmetric bilayer results compatible with the anomalous 
C* extracted for the incomplete bilayer. 

In his study, Troyer fitted a line to QMC data in the 
range 0.1 < T/J < 0.4 [ll|, resulting in the dashed 
line reproduced in Fig. [3] — the T = intercept is con- 
sistent with 1/4 within statistical errors. However, the 
line deviates considerably from our size-converged results 
for T/J > 0.1. The discrepancy may be due to a dif- 
ferent size- convergence procedure; an 1 / L extrapolation 
was mentioned in Ref. [ll| . whereas we have used a cri- 
terion of size-independence. The latter procedure should 
be more reliable because an asymptotically exponential 
convergence is expected at finite temperature. An l/L 
extrapolation can lead to a too high value when fitting 
only to a range of points for which an almost linear in 
l/L behavior is observed, as can be seen in the inset of 
Fig. 1 for T/J = 0.25. 

Wc now turn to our results for an S — 1 impurity, 
which we have realized in the symmetric bilayer as shown 
in Fig. [BJc). As shown in Fig. the finite-size behavior 
of the impurity susceptibility is similar to the S — 1/2 
impurities. For L = 4, 8, and 16 the asymptotic low T 
behavior T X f mp (T -4 0) = S(S + l)/3 = 2/3 is clearly 
seen. Size converged low-T results are shown in the in- 
set. A linear behavior sets in at T/J as 0.15. Here the 
extrapolated C* — 0.663(2) is slightly below the normal 
value 2/3, but the difference is too small to definitely con- 
clude that this is the case. The results do show that the 
anomaly is smaller than for S = 1/2. The theory does 



not predict how the fractional spin evolves as a function 
of 5 Si- 

In summary, we have presented evidence from unbiased 
numerical computations of an anomalous Curie response 
of an S = 1/2 impurity spin in a 2D quantum-critical 
antiferromagnet. The anomalous Curie constant C* = 
0.262(2) is only w 5% larger than the normal C = 1/4. 
For S — 1 we obtain C* marginally below the normal 
value, 2/3, but better statistics is needed to confirm this. 

It should be noted that the anomalous Curie constant 
cannot be strictly interpreted as due to a fractionalized 
impurity spin, because it is a finite-temperature quantity 
with contributions from many states with different total 
spin, even as T — > 0. Fractionalization does not occur in 
the ground state for finite L, but an interesting universal 
impurity-induced spatial structure has been found [17 1. 
Recently impurity effects have been examined theoreti- 
cally also in fractionalized spin liquid states ID, 19j |. 
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